library(magrittr)
It should be straightforward to use DBSLMM outputs and fam files to write custom code for CV+.
Inputs needed:
alpha <- 0.1
out <- list()
res_file <- "res-continuous.rds"
if (!file.exists(res_file)){
res <- lapply(1:25, function(x) {
knitr::knit_child(
'cv-plus-child.Rmd',
envir = environment(),
quiet = TRUE
)
get("out")
})
saveRDS(res, file = res_file)
} else {
res <- readRDS(res_file)
}
coverages <- lapply(X = res,
FUN = function(mytib){
mytib %>%
dplyr::summarise(coverage = mean(in_interval))
}
)
coverage_tib <- do.call("rbind", coverages) %>%
tibble::as_tibble() %>%
dplyr::mutate(trait_num = 1:length(res))
coverage_tib %>%
gt::gt()
| coverage | trait_num |
|---|---|
| 0.774 | 1 |
| 0.802 | 2 |
| 0.832 | 3 |
| 0.826 | 4 |
| 0.868 | 5 |
| 0.826 | 6 |
| 0.878 | 7 |
| 0.834 | 8 |
| 0.842 | 9 |
| 0.834 | 10 |
| 0.864 | 11 |
| 0.886 | 12 |
| 0.874 | 13 |
| 0.886 | 14 |
| 0.858 | 15 |
| 0.846 | 16 |
| 0.884 | 17 |
| 0.878 | 18 |
| 0.838 | 19 |
| 0.888 | 20 |
| 0.910 | 21 |
| 0.838 | 22 |
| 0.820 | 23 |
| 0.838 | 24 |
| 0.804 | 25 |
for (i in 1:length(res)){
res[[i]] <- res[[i]] %>%
dplyr::mutate(trait_num = i)
}
dplyr::bind_rows(res) %>%
ggplot2::ggplot() +
ggplot2::geom_boxplot(ggplot2::aes(x = trait_num, y = interval_width, group = trait_num))
h_tib <- readr::read_csv("~/research/ukb-intervals/shell_scripts/h2.csv", col_names = FALSE) %>% dplyr::select(-26) # drop 26, which is due to presence of trailing commas
## Rows: 5 Columns: 26
## ── Column specification ───────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (25): X1, X2, X3, X4, X5, X6, X7, X8, X9, X10, X11, X12, X13, X14, X15, ...
## lgl (1): X26
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
pp <- h_tib %>%
dplyr::summarise(dplyr::across(X1:X25, mean)) %>%
tidyr::pivot_longer(X1:X25) %>%
dplyr::rename(mean_herit = value) %>%
dplyr::mutate(trait_num = 1:25) %>%
dplyr::select(-name) %>%
dplyr::select(trait_num, mean_herit) %>%
dplyr::left_join(dplyr::bind_rows(res), by = "trait_num")
pp2 <- pp %>%
ggplot2::ggplot() + ggplot2::geom_boxplot(ggplot2::aes(x = mean_herit, y = interval_width, group = trait_num, colour = as.factor(trait_num), label = trait_num), show.legend = FALSE)
## Warning: Ignoring unknown aesthetics: label
pp2 %>% plotly::ggplotly(tooltip = "label")
pp2 <- pp %>%
dplyr::group_by(trait_num) %>%
dplyr::summarise(mean_herit = mean(mean_herit), coverage = mean(in_interval)) %>%
ggplot2::ggplot() + ggplot2::geom_point(ggplot2::aes(x = mean_herit, y = coverage, colour = as.factor(trait_num)))
pp2 %>%
plotly::ggplotly()
summ <- pp %>%
dplyr::group_by(trait_num) %>%
dplyr::summarise(mean_herit = mean(mean_herit), coverage = mean(in_interval))
cor(summ$mean_herit, summ$coverage)
## [1] -0.6851291